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This  report  presents  a  formal  approach  for  optimizing  the  shape  of  freely  supported  isotropic 
plates  to  withstand  air  blast  loading.  Unique  difficulties  are  presented  by  the  nature  of  short- 
duration  dynamic  loading,  viz.  transient  dynamic  response,  monitoring  of  maximum  plastic 
strain  failure  at  every  point  in  the  panel  over  time,  optimizers  that  can  handle  non-differentiable, 
nonconvex  and  computationally  expensive  functions,  and  mesh  distortion.  LS-DYNA  is  used  as 
the  finite  element  software.  The  finite  element  model  has  been  developed  to  reflect  experimental 
test  conditions  and  observed  structural  response.  The  goal  has  been  to  minimize  dynamic 
displacement  relative  to  the  fixture,  while  monitoring  plastic  strain  values,  mass,  and  envelope 
constraints.  A  Fortran  code  has  been  developed  to  implement  the  methodology.  Sinusoidal  basis 
shapes  are  used  to  obtain  an  optimized  double-bulge  shape.  Importantly,  the  flat  plate  is 
associated  with  a  concentration  of  plastic  strain  at  its  center  while  for  the  optimized  bulge  shape, 
the  plastic  strain  is  smeared  around  the  support  showing  greater  utilization  of  material.  Results 
show  that  much  superior  structural  systems  can  be  designed  compared  to  ad-hoc  techniques  that 
sometimes  fail  to  improve  even  the  baseline  design  of  a  flat  plate.  Change  in  optimized  shape 
with  increasing  offset  in  charge  location  is  studied.  A  methodology  for  optimizing  honeycomb 
sandwich  structures  is  presented,  utilizing  a  novel  technique  for  linking  honeycomb  core 
geometry  with  its  stress-strain  curve. 
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Nomenclature 

x  =  design  variables  vector 

xL  =  lower  limit  on  design  variables 

xu  =  upper  limit  on  design  variables 

G  =  vector  of  x-,  y-,  z-  coordinates  of  nodes  in  FE  model 

q'  =  i  th  velocity  field  or  trial  shape  change  vector 

w  =  z-  or  normal  displacement 

£max  =  maximum  plastic  strain  at  failure  for  the  material 

s  =  equivalent  plastic  strain  vector 

M  =  total  mass  of  the  structure 

Mm ax  =  upper  limit  for  the  mass  of  the  structure 

t  =  thickness  of  the  structure  (plate)  at  any  (pc,  y)  location  in  the  plate 

tmm  =  Minimum  thickness  allowed 

z  =  vector  of  z-coordinate  of  the  nodes 

zu  =  upper  limit  on  z-coordinate 

zL  =  lower  limit  on  z-coordinate 

(det  Jj)  =  Jacobian  of  j  th  hexahedral  element  at  all  the  eight  nodes 
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1.  Introduction 

This  report  presents  an  approach  and  methodology  to  optimize  the  shape  of  panels  to 
mitigate  the  effect  of  air  blast  loading,  which  comprises  of  a  short  duration  pressure  pulse. 
Unique  difficulties  are  presented  by  the  nature  of  the  blast  loading,  viz.  transient  dynamic 
response,  modeling  the  experimental  setup  with  appropriate  boundary  conditions,  monitoring  of 
maximum  plastic  strain  failure  at  every  point  in  the  panel  over  time,  and  optimizers  that  can 
handle  non-differentiable,  nonconvex,  and  computationally  expensive  functions.  Much  superior 
structural  shapes  can  be  designed  using  a  formal  optimization  methodology  compared  to  using 
ad-hoc  techniques. 

Considerable  attention  in  journals  and  conferences  worldwide  has  been  given  to  analysis 
of  metallic  and  composite  panels,  subject  to  both  blast  and  ballistic  loads.  Regarding  designing 
for  impact  mitigation,  much  greater  focus  has  been  placed  on  ballistic  impact  rather  than  on 
blast.  Very  few  papers  use  formal  optimization  techniques.  Publications  relating  to  blast  damage 
mitigation  are  given  below. 

This  report  focuses  on  shape  of  isotropic  metal  plates.  Dharaneepathy  and  Sudhesh  [1] 
investigate  stiffener  patterns  on  a  square  plate  subject  to  blast  loads  modeled  using  Friedlander’s 
exponential  function.  While  formal  optimization  was  not  used  to  optimize  the  stiffener  patterns, 
they  demonstrated  that  stiffeners  do  provide  significant  advantage  compared  to  an  unstiffened 
panel  of  same  weight,  and  that  a  waffle  pattern  is  not  as  good  as  a  novel  pattern  that  they 
proposed.  Failure  was  not  considered  in  their  study  -  that  is,  only  deflection  was  considered.  Hou 
et  al  [2]  also  consider  the  effect  of  stiffener  size  on  blast  response.  Xue  and  Hutchinson  [3]  and 
Fleck  and  Deshpande  [4]  compare  blast  resistance  of  solid  versus  sandwich  panels  (such  as 
pyramidal  truss  core,  square  honeycomb  and  folded  plate).  The  plates  were  considered  to  be 
infinitely  long  in  one  direction  and  fixed  at  the  ends  of  the  short  direction.  ABAQUS/Explicit 
was  used  to  model  the  blast  load  in  Ref.  [3]  while  an  approximate  analytical  approach  was  used 
in  Ref.  [4],  Blast  in  both  air  and  water  were  included  in  the  comparative  study.  From  certain 
normalized  displacement  versus  impulse  graphs,  it  was  concluded  that  some  of  the  sandwich 
topologies  outperformed  solid  panels  of  same  mass,  especially  in  water.  Yen,  Skaags  and 
Cheeseman  [5]  present  an  experimentally  validated  dynamic  analysis  procedure  utilizing  Ls- 
Dyna  and  the  ConWep  air  blast  function  with  shock  mitigation  materials  such  as  honeycomb  or 
foam.  The  numerical  results  indicate  that  significant  reduction  in  the  maximum  stress  amplitude 
propagating  within  the  protected  components  can  be  achieved  by  suitable  selection  of 
honeycomb  material  with  proper  crush  strength.  Liang,  Yang  and  Wu  [6]  focused  on  static 
loading  and  hence  not  very  relevant  in  the  present  context  of  blast  effects  where  inertia  effects 
play  a  role.  Further,  their  paper  ignores  material  nonlinearity.  Main  and  Gazonas  [7]  study  the 
effect  of  an  air  blast  on  uniaxial  crushing  of  a  cellular  sandwich  plates.  They  investigated  the 
physics  of  shock  mitigation  for  different  geometries  of  the  cellular  core.  Icardi  and  Ferrero  [8] 
study  optimum  fiber  orientations  in  a  laminated  composite  to  absorb  energy  while  maintaining 
stiffness.  Further  details  on  effectiveness  of  blast  mitigation  solutions  in  a  laboratory  using  the 
Pendulum  Test  are  given  in  [5,  9,  10]. 

There  are  instances  in  the  public  domain  that  show  the  mitigating  effects  of  the  shape  of 
the  plate.  For  example,  United  States  Patent  7357062 

(http://www.freepatentsonline.eom/7 357 062 .html )  shows  that  V-shape  deflects  blast  waves 
away.  We  view  the  V-shape  as  a  bulge  towards  the  charge  and  belonging  to  the  family  of  single- 
and  double-  bulges  obtained  in  this  formal  optimization  study.  Our  focus  is  on  the  development 
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of  a  methodology  that  can  be  applied  to  the  design  of  any  plate,  i.e.  not  just  aluminum  plates,  and 
design  for  any  given  amount  of  blast  charge,  location  and  any  boundary  conditions. 

In  contrast  with  the  papers  cited  above,  this  report  uses  formal  optimization  techniques, 
focuses  on  solid  metal  plates,  uses  Ls-Dyna  for  simulation  and  considers  plastic  strain  to  failure 
in  the  metal.  The  methodology  involves  integrating  an  optimizer  with  the  Ls-Dyna  simulation 
code.  In  Section  2,  we  formulate  the  optimization  problem.  Ls-Dyna  finite  element  modeling  of 
blast  loading  is  discussed  in  Section  3.  Material  properties  are  also  given  in  Section  3.  Sections 
4-5  contain  the  shape  optimization  approach.  Computer  code  development  is  discussed  in 
Section  6,  results  in  Section  7  and  conclusions  in  Section  8.  Clear  improvement  over  a  flat  plate 
is  demonstrated. 

2.  Problem  Definition 

The  schematic  diagram  of  the  plate  used  for  shape  optimization  is  shown  in  Fig.  1.  The 
standoff  distance  of  the  charge  is  taken  to  be  0.4064  m.  It  should  be  noted  that  plate  is  just  a  part 
of  a  freely  supported  ‘grip’  assembly  used  to  model  the  experimental  condition  as  explained  in 
the  next  Section.  With  reference  to  Fig.  1,  the  basic  problem  addressed  in  this  work  can  be  stated 
as  follows: 

Given  a  set  of  basis  shapes  which  controls  the  shape  and  thickness  of  the  plate,  a  mass  limit  for 
the  structure,  plastic  strain  limits  representing  fracture  strength,  a  minimum  thickness  for  the 
panel,  and  a  geometric  envelope  within  which  the  structure  must  lie ,  determine  the  best  possible 
combination  of  these  basis  shapes  that  minimizes  the  deflection  (at  the  first  peak  in  time). 


Figure  1.  Baseline  Structural  Model:  Plate  dimensions  are 

(1.2192  m  x  1.2192  m  x  0.024  m),  Standoff  distance  d  =  0.4064  m 
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Utilizing  the  already  described  notation,  the  problem  can  be  formulated  as  follows. 


minimize 

Hx)|, 

(1) 

subject  to 

£/'  —  '-'max 

M  <  Mmax 

for  each  element  j 

t  —  t  m  j  n 

xL<x<xu 
det  Jj  (x)  >  0 

for  each  element  j 

zL  <  z  <  zu 

(envelope) 

Important  notes  regarding  the  above  problem  are  as  follows. 

(i)  Euclidean  norm  of  the  relative  z-displacement  of  nodes  in  the  plate  is  taken  as  the 
objective  function  (i.e.  r  =  2  in  the  objective  function  definition  above),  x-  and  y- 
displacements  are  not  significant  and  are  not  considered.  The  term  ‘relative’  is 
explained  subsequently.  The  displacement  is  a  function  of  time,  and  value  at  first 
peak  is  monitored. 

(ii)  Plastic  strain,  also  a  function  of  time,  stabilizes  after  a  certain  simulation  time 
duration.  This  stabilized  value  is  used  in  the  constraint. 

(iii)  M  above  refers  to  the  combined  mass  of  the  assembly  (Fig.  2). 

(iv)  Thickness  is  computed  from  nodal  coordinates  of  the  hexahedral  elements  used  in 
the  FE  model.  Element  distortion  is  prevented  by  computing  determinant  of 
Jacobian  in  every  element  and  is  forced  to  stay  positive  during  optimization.  The 
Jacobian  for  each  element  is  computed  from  element  nodal  coordinates. 

(v)  The  focus  is  on  a  freely  supported  plate  -  the  methodology  can  be  readily  applied 
to  plates  with  other  boundary  conditions  noting  that  the  optimum  shapes  may  be 
different. 

3.  Ls-Dyna  Modeling  Considerations 

Initially,  a  uniformly  thick  (i.e.  flat)  square  aluminum  plate  is  considered  as  a  starting 
shape  or  baseline  design.  Initial  studies  were  carried  out  with  the  baseline  design  to  better 
understand  the  explicit  finite  element  analysis  procedure. 

Boundary  conditions  study:  Various  types  of  boundary  conditions  along  the  plate  edges  were 
studied  for  a  freely  supported  plate,  including  one  or  more  ‘fences’  of  springs  around  the  edges. 
While  the  plate  on  springs  is  discussed  in  Section  7,  all  these  modeling  attempts  produced 
unrealistic  plastic  strain  concentrations  at  the  comers.  A  grip  system  or  assembly  that  holds  the 
plate  was  finally  adopted  owing  to  best  correlation  with  experimentally  observed  plate  failures 

(Fig-  2). 
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Rigid  plates 


Plate  of 
interest 

Aluminum 
tiller  plate 


Figure  2.  Exploded  view  of  the  structure  modeled 

The  whole  assembly  is  free  to  move  during  blast.  That  is,  the  assembly  is  unconstrained,  to 
closely  simulate  blast  testing  in  the  field  and  in  pendulum  based  experiments  as  noted  above. 
Appropriate  boundary  conditions  using  contact  surfaces  are  used  to  model  the  assembly.  The 
simulation  on  the  baseline  (i.e.  flat  plate)  design  produces  maximum  plastic  strain  at  the  center  of 
the  plate  as  also  observed  in  testing.  Referring  to  Fig.  2,  components  of  the  grip  assembly  are: 

i)  the  ‘plate  of  interest’  (red),  made  of  aluminum. 

ii)  a  filler  plate  (blue),  also  made  of  aluminum,  having  same  thickness  as  the  edges 
of  the  plate  of  interest,  necessary  for  the  assembly. 

iii)  three  rigid  plates,  whose  material  properties  are  not  relevant  -  steel  is  assumed  for 
density,  and  an  artificially  very  high  value  for  Young’s  modulus  is  used  to  ensure 
rigidity.  Importantly,  the  circular  cover  plate  (blue)  above  the  plate  of  interest 
creates  peak  plastic  strain  at  the  center  of  the  flat  plate  of  interest,  as  this  is  the 
location  where  a  freely  supported  flat  plate  is  experimentally  observed  to  fail. 
That  is,  without  this  circular  cover  plate,  as  with  other  modeling  attempts 
mentioned  above  such  as  spring  supports,  unrealistic  pseudo  concentrations  of 
plastic  strains  are  observed  at  comers. 

iv)  A  square  portion  of  the  model  at  the  center  of  plate  of  interest  (red)  is  taken  to  be 
the  domain  of  shape  optimization.  This  square  portion  lies  within  the  circle  of  the 
cover  plate  (blue).  Thus,  shape  changes  (such  as  bulges)  are  introduced  in  this 
square  domain  only.  The  region  outside  the  domain  in  the  plate  of  interest  is  flat, 
whose  thickness  equals  that  of  the  surrounding  (blue)  filler  plate. 

It  is  important  to  get  boundary  conditions  to  model  reality  as  the  optimum  shape  will  depend  on 
it.  To  eliminate  the  rigid  body  component,  the  w(x)  in  the  objective  function  in  Eq.  (1)  is  taken 
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to  be  the  ‘relative’  displacement  that  is  obtained  by  subtracting  the  nodal  displacement  at  a  point 
in  the  plate  from  the  displacement  of  a  reference  point  in  the  rigid  grip. 

Elements:  All  components  are  modeled  using  8-noded  hexahedral  or  brick  elements.  The  ‘plate’ 
is  thus  more  a  ‘solid’  since  brick  elements  and  not  shell  elements  are  used. 

Mesh  sensitivity  study:  Further,  it  is  desirable  to  establish  a  surrogate  model  that  is 
computationally  efficient  and  is  accurate  enough  to  capture  the  response  trends  correctly  for 
iterative  optimization.  With  respect  to  this  problem,  the  displacement  and  plastic  strain 
predictions  are  important.  Different  mesh  densities  were  chosen  in  this  study.  A  study  was 
conducted  with  different  mesh  densities  ranging  from  4862  elements  to  51902  elements.  It  was 
found  that  peak  displacement  was  not  sensitive  to  mesh  refinement,  while  maximum  plastic 
strain  was.  The  FE  analysis  time  for  one  complete  analysis  varied  between  90  s  for  the  4862- 
element  model  to  15  min  for  the  51902-element  model  when  runs  were  made  on  a  Intel  P4-3.6 
GHz  machine  with  3  GB  RAM.  The  simulation  time  is  taken  to  be  the  time  required  for  plastic 
strain  to  stabilize  at  a  constant  value.  The  coarser  FE  model  was  chosen  here  with  appropriately 
reduced  plastic  strain  limit  for  constraint  evaluation  in  the  optimization  problem.  Upon  solution, 
a  fine  mesh  of  the  optimized  shape  was  used  to  check  that  plastic  strains  were  less  than  actual 
failure  limits. 

Loading:  Blast  load  exerted  is  calculated  using  the  ConWep  function  in  Ls-Dyna  The  inputs  for 
*LOAD_BLAST  command  are  equivalent  TNT  mass,  type  of  blast  (surface  or  air),  load  curve, 
charge  location,  and  surface  identification  for  which  pressure  is  applied.  ConWep  calculates  the 
appropriate  reflected  pressure  values  and  then  applies  these  to  the  appropriate  surfaces  by  taking 
account  the  angle  of  incidence  of  the  blast  wave.  It  should  be  noted  that  the  loading  changes 
with  changes  in  shape  of  the  plate,  although  this  is  automatically  calculated  during  each  finite 
element  analysis.  Blast  parameters  are  given  in  Table  1. 


Table  1.  Blast  Load  Input  Data 


Property 

Value 

Equivalent  mass  of  TNT 

1  kg 

Blast  Location 

(0.0,0.0,-0.4064)  m 

Type  of  Burst 

Air  Blast  (Spherical 

Charge) 

Material  model:  The  to-be-designed  plate  and  the  surrounding  filler  plate  are  made  of 
Aluminum  5083  with  the  bilinear  elastic-plastic  material  model.  Material  properties  are  listed  in 
the  Table  2.  The  *MAT_PLASTIC_KINEMATIC  input  card  is  used. 


7 


Table  2.  Aluminum  5083  Material  Properties 


Property 

Value 

Mass  Density 

3 

2700  kg/m 

Young’s  Modulus 

68.9  GPa 

Poisson’s  Ratio 

0.33 

Yield  Stress 

225  GPa 

Tangent  Modulus 

633  GPa 

Hardening  Parameter 

1.0 

Failure  Strain 

0.39 

4.  Shape  Optimization  Methodology 

Prior  to  using  formal  optimization,  ad-hoc  studies  were  carried  out  to  generate  improved 
shapes.  The  following  designs  were  studied: 

a)  a  singly  corrugated  plate,  and 

b)  a  dimpled  (indented)  plate,  with  dimples  facing  towards  and  away  from  the  charge, 
respectively. 

These  ad-hoc  attempts,  however,  did  not  produce  any  reduction  in  maximum  deflection  and, 
moreover,  lead  to  rupture  levels  of  plastic  strains.  Hence,  it  was  decided  to  use  formal 
optimization  methodology. 

A  square  portion  of  the  model  at  the  center  of  plate  is  taken  to  be  the  domain  of  shape 
optimization.  That  is,  shapes  changes  only  occur  in  this  region.  This  ensures  that  changes  in  plate 
shape  do  not  result  in  changes  in  the  grip  system  (see  the  circular  blue  cover  plate  in  Fig.  2). 
However,  a  thickness  change  in  the  design  plate  has  to  be  matched  by  an  equal  thickness  change 
in  the  filler  plate  for  the  assembly  to  function. 

The  key  equation  to  implement  shape  optimization  is  [1 1,  12] 

G(x)  =  Goriginal  +  q'  (2) 

/=i 

where  G  is  a  grid  point  coordinates  vector,  representing  x-,  y-,  z-  coordinates  of  all  nodes  in  the 
model.  Each  x/f  represents  the  amplitude  of  a  ‘permissible  shape  change  vector’  or  what  is 
commonly  called  a  ‘velocity  field’  vector  q\  Velocity  fields  have  nothing  to  do  with  actual 
velocities  of  the  model  under  loading.  Vectors  {q*}  are  generated  outside  the  iterative 
optimization  loop.  Goriginai  is  the  current  (flat)  shape.  Visualization  of  a  {q'j  is  identical  to 
visualization  of  a  displacement  field  in  finite  elements:  {q4}  is  multiplied  by  a  magnification 
scalar  and  added  to  the  current  grid  to  obtain  a  displaced  grid,  except  that  here  the  displaced  grid 
represents  a  new  shape  and  is  called  a  basis  shape. 

The  role  of  the  optimizer  is  to  choose  x  so  that  the  corresponding  shape  G(x  )  is 
optimum.  As  x  is  iteratively  changed  by  the  optimizer,  the  grid  point  coordinates  G  are  updated, 
a  FE  input  file  is  then  written  and  an  analysis  is  carried  out  to  evaluate  the  various  functions  in 
the  optimization  problem.  This  flow  of  information  may  be  shown  as  follows. 
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Optimizer  calls  with  a  xk  at  the  Ath  iteration  ->  construct  shape  G  from  Eq.  (2)  ->  Ls-Dyna 
analysis  ->  evaluate  objective  and  constraint  functions  ->  return  to  optimizer 


Chosen  Basis  Shapes  (or  Velocity  Fields) 

Sinusoidal  velocity  fields  for  the  top  and  bottom  surfaces,  independently,  are  chosen  based  on 
the  following  equation: 

.  _  .  m  7i x  .  rury 

,72)=  C  sin -  sm -  (3) 

L  L 

where  C  is  a  suitable  normalization  factor  and  m,  n  are  integers  taking  on  values  1,2,3,  etc.  In 
this  report,  optimization  results  are  presented  for  two  cases  -  three  design  variables  (denoted  3- 
DV  Case)  and  nine  design  variables  (denoted  9-DV  Case). 

3  -DV  Case :  m  =  n  =  1.  This  gives  a  total  of  three  (3)  symmetric  basis  shapes  corresponding  to 
q1  =/(l,l)  for  the  top  surface, 
q2  =/(l,l)  for  the  bottom  surface, 
q3  =  thickness  change 

Specifically,  q1  represents  a  bulge  in  the  shape  of  the  top  surface  while  bottom  surface  is  fixed 
(other  thru-thickness  nodes  are  moved  to  preserve  equal  spacing),  q2  represents  a  bulge  on  the 
bottom  surface  while  top  surface  is  fixed,  and  q3  =  a  thickness  change  only  (that  is,  middle  layer 
of  nodes  in  the  plate  are  fixed  while  top  and  bottom  surfaces  move  in  opposite  directions).  Thus, 
the  design  variable  vector  is  x  =  [xi  ,  X2,  *3  ]T  ,  and  the  optimizer  tries  to  determine  an  optimum 
combination  of  these  three  basis  shapes.  Basis  shape  corresponding  to  q1  is  illustrated  in  Fig.  3. 
Note  that  the  bulge  can  be  positive  or  negative  depending  on  the  sign  of  xt 


(a)  (b) 


Figure  3.  (a)  Basis  shape  corresponding  to  top  surface  bulge,  (b)  Basis  shape  as  part  of  the 
full  plate,  showing  the  domain  for  shape  optimization 

2 -DV  Case:  a  special  case  of  the  3-DV  case  where  only  q1  and  q2  are  used.  That  is,  there  is  no 
change  in  the  plate  edge-thickness.  However,  the  thickness  in  the  interior  of  the  plate  will  be 
non-uniform  as  these  are  determined  by  the  relative  magnitudes  of  the  bulges  created  by  q1  and 
q2.  Note  that  the  plate  is  modeled  using  3-D  hexahedral  elements  and  not  plate  elements. 

9-DV  Case:  m  =  n  =  2.  This  gives  a  total  of  nine  (9)  basis  shapes  corresponding  to 
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q1  =/(  1,1),  q2  =/(l,2),  q3  =f( 2,1),  q4  =f{2,2)  for  the  top  surface, 
q5  =/(  1,1),  q6  =/(l,2),  q7  =/(2,l),  q8  =f(2,2)  for  the  bottom  surface, 
q9  =  thickness  change 

Thus,  x  is  a  (9x1)  vector.  It  should  be  notes  that  unlike  the  3-DV  Case  above,  presence  of 
unsymmetric  basis  shapes  may  lead  to  an  unsymmetric  final  optimum  shape. 

5.  Optimizers 

During  computational  experiments,  it  was  observed  when  using  a  gradient-based 
optimizer,  downhill  or  descent  search  directions  did  not  always  lead  to  a  reduction  in  the 
objective  function  even  for  small  steps.  The  problem  was  seen  to  be  clearly  due  to  non- 
differentiable  functions,  attributable  to  the  dynamic  nature  of  the  response.  Hence  the  use  of 
gradient  based  optimizers  is  not  appropriate.  The  Differential  Evolution  (DE)  technique  has 
proven  to  be  successful  [13],  DE  is  similar  to  genetic  algorithms  in  some  respects  such  as 
involving  a  population  of  designs  and  having  generations  of  designs.  DE  requires  fewer  control 
variables,  is  robust  and  is  very  well  designed  for  parallel  computation  implementation.  Decision 
parameters  of  the  algorithm  are  mutation  scaling  factor  and  cross  over  factor  for  the  generation 
of  a  population  during  a  new  generation.  Here,  random  scaling  factor  is  used  for  the  linear 
crossover  combination  of  best  member  and  older  population,  for  better  diversity.  Penalty 
approach  is  used  to  satisfy  the  constraints.  Quadratic  penalty  is  used  for  plastic  strain,  geometric 
and  mass  constraints,  while  violation  of  Jacobian  constraint  is  handled  by  returning  a  very  high 
function  value  since  finite  element  analysis  cannot  be  carried  out  with  a  distorted  mesh.  The  total 
number  of  finite  element  analysis  is  the  product  of  population  size  and  the  number  of 
generations.  Since  DE  is  stochastic  in  nature,  different  seeds  have  been  tried  and  the  best  answer 
based  on  the  minimum  objective  function  value  is  chosen. 

Other  optimizers  such  as  a  genetic  algorithm  have  also  been  used.  But  the  DE  code 
performed  better  for  this  particular  problem.  The  use  of  design  of  experiments  (DOE)  combined 
with  response  surface  optimization  may  be  another  route  to  take,  however,  the  nonlinearity  of 
the  problem  will  entail  several  repeated  best- fit  surfaces  followed  by  optimization  within  some 
iterative  scheme.  In  any  vase,  parallelization  of  the  DE  optimizer  makes  it  viable  to  directly 
optimize  with  accurate  Ls-Dyna  FE  model  response. 


Parallel  Computations 

We  set  the  population  size  equal  to  ten  times  the  number  of  design  variables.  Thus,  in  the 
9-DV  case,  population  size  is  90.  With  the  number  of  generations  set  as  40,  we  have  a  total  of 
3600  finite  element  analyses  or  Ls-Dyna  runs.  To  reduce  the  total  execution  time,  a  simple 
parallelization  of  the  code  using  MPI  calls  has  been  implemented  [14],  We  used  the  LION-XC 
cluster  at  Penn  State  University.  Reduction  in  the  total  time  compared  to  execution  on  a  single 
processor  workstation  makes  this  optimization  approach  feasible  in  the  9-DV  case.  LION-XC  is 
a  cluster  with  each  compute  node  being  a  dual  3.0-GHz  Intel  Xeon  3160  (Woodcrest)  Dual-Core 
Processors 
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6.  Computer  Code  Development 

A  Fortran  computer  program  has  been  developed  to  implement  the  above  procedures. 
The  main  blocks  are  shown  in  Fig.  4. 


Figure  4.  Flow  Diagram  of  the  computer  code  for  shape  optimization  using  Ls-Dyna 


A  key  element  of  the  implementation  lies  in  the  creation  of  two  types  of  input  files,  as  discussed 
below. 

(A)  Input  File  for  Design  Optimization  contains: 

(i)  seed  used  for  DE, 

(ii)  generation  limit, 

(iii)  population  size, 

(iv)  geometric  envelope  limits  (limits  for  the  z-coordinate  of  the  nodes), 

(v)  maximum  mass  limit, 

(vi)  plastic  strain  limit, 

(vii)  upper  and  lower  limit  of  design  variable, 

(viii)  velocity  fields,  and 

(ix)  model  related  inputs. 

(B)  Input  File  for  FE  Analysis  (Ls-Dyna)  contains: 

(i)  charge  data, 

(ii)  material  properties,  and 

(iii)  nodal  coordinates  and  element  connectivity. 

Typical  values  used  are  given  in  the  Table  3.  Note  that  plastic  strain  limit  is  less  than  the  material 
property,  to  account  for  the  fact  that  a  coarser  mesh  is  being  used  for  optimization  (based  on  the 
discussion  in  Section  3).  Coordinate  files  are  written  for  every  population  generated  and  checked 
for  distortion  of  the  mesh.  High  function  value  is  returned  for  a  distorted  mesh,  otherwise  the  Ls- 
Dyna  solver  is  invoked.  Ls-Dyna  writes  nodal  and  element  related  time  history  outputs  to  ASCII 
files  called  ‘nodout’  and  ‘elout’  respectively.  Fortran  routines  are  written  to  open  these  files  and 
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calculate  relative  displacement  values  and  constraint  values  which  are  added  to  objective 
function  as  quadratic  penalties.  Based  on  the  objective  function  values  of  all  the  population 
members  and  previous  generations,  best  member  is  selected  and  stored.  After  the  last  generation, 
best  member  is  written  to  the  output  file  along  with  max  displacement  value  and  values  of  plastic 
strain.  The  visualization  of  the  results  is  through  LS-Prepost. 

Table  3.  Typical  values  of  input  parameters  used  in  the  input  file 


Parameter 

Value 

Generation  limit 

20  to  45 

Population  size 

30  to  200 

Envelope  limits  ( m ) 

-0.3  to  1.0  (large) 

-0.1  to  0.14  (small) 

Max  mass  of  assembly 

(kg) 

1890.0 

Plastic  strain  limit 

0.15 

Scaling  factor,  C  in  Eq.  (2) 

0.02 

Seed 

1170,  2349 

7.  Results  -  Optimized  Panel  Shapes 

Results  using  2-DV,  3-DV  and  9-DV  (described  in  Section  4)  are  presented  below 
followed  by  physical  interpretations  and  some  sensitivity  studies.  Small  variations  in  the 
optimized  mass  from  the  limiting  value  is  owing  to  use  of  exterior  penalty  functions  in  the 
stochastic  optimizer. 

2-DV  Case:  Results  are  presented  in  Table  4,  where  the  edge -thickness  of  the  plate  is  held 
constant,  using  only  two  velocity  fields  q1,  q2.  Since  edge-thickness  is  constant,  the  surrounding 
filler  plate  thickness  is  also  unchanged  and  does  not  enter  into  the  problem  during  the  iterative 
optimization  process.  A  single-bulge  towards  the  charge  is  obtained  as  optimum  (Fig.  5). 

Comparing  initial  and  optimized  designs,  RMS  displacement  is  reduced  from  20.5 1  mm 
to  5.84  mm,  corresponding  to  a  72%  improvement,  ft  may  be  argued  that  this  comparison  is  not 
accurate  as  there  is  a  difference  due  to  use  of  penalty  functions  in  the  optimizer  initial  (baseline) 
mass  and  optimized  mass,  viz.  between  1872.2  kg  and  1899.6  kg  (a  1.5%  difference).  Thus,  for  a 
more  accurate  comparison,  the  optimized  plate  is  also  compared  to  a  flat  plate  of  equal  mass. 
This  has  been  done  by  adjusting  the  plate  thickness.  The  total  mass  of  the  structure,  is  maintained 
at  1899.5  kg,  to  avoid  changes  to  total  baseline  impulse,  by  adjustment  of  the  density  of  the  rigid 
components  whose  properties  are  not  relevant  except  with  regard  to  providing  a  rigid  grip  to  the 
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plate.  It  is  found  that  RMS  displacement  is  reduced  by  67%  (compared  to  72%  quoted  for 
baseline),  plastic  strain  is  reduced  by  80%  (instead  of  84%)  and  impulse  is  reduced  by  13.3% 
(compared  to  13.7%).  Results  for  both  a  large  envelope  (LE)  and  a  small  envelope  (SE)  are 
presented.  The  Targe  envelope’  case  is  slightly  better  compared  to  ‘small  envelope’  case,  as  is  to 
be  expected.  Note  that  the  envelope  is  a  user-defined  constraint. 


Table  4.  Optimum  solution  for  2-DV  with  Small  and  Large  Envelope 


Property 

Initial  Flat 

(Baseline) 

Plate 

Optimized 

2-DV  LE 

Optimized 

2-DV  SE 

Population  Size  in  DE  optimizer 

30 

30 

Number  of  Generations  in  DE 

45 

45 

Max.  relative  displacement  (mm) 

occurs  at  1st  peak 

58.43 

10.5 

14.28 

Objective  =  RMS  displacement,  mm 

20.51 

5.84 

7.34 

Max.  plastic  strain 

0.1277 

0.021 

0.021 

Total  mass  of  structure  (kg) 

1872.2 

1899.6 

1903.81 

Total  Z-momentum  (kN-sec) 

(i.e.  saturated  impulse) 

6.254 

5.40 

5.41 

2-DV  SE:  Side  View 

2-DV  SE:  Isometric  View 

2-DV  LE:  Side  View 

2-DV  LE:  Isometric  View 

Figure  5.  Final  optimum  shapes  with  2-DV  where  plate  edge  thickness  unchanged 

(a  single  bulge  towards  the  charge) 


13 


3-DV  Case:  Three  basis  shapes  are  used  here.  Initial  baseline  (flat  plate)  model  and  the  resulting 
optimum  shape  of  the  panel  are  shown  in  Figs.  6a  and  6b.  The  optimized  shape  turns  out  to  be  a 
double  bulge  (Figs.  6b,  7)  rather  than  the  single  bulge  in  the  2-DV  case. 


<a>  (b) 

Figure  6.  (a)  Baseline  design(uniform  thickness) ,  (b)  Optimized  plate  for  large 
envelope  case:  double  bulge 


1 

0.0844m 
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Figure  7.  Dimension  of  the  optimal  design  panel  (3-DV  large  envelope  case) 


Thickness  of  the  plate  of  interest  outside  the  shape  domain  region,  which  is  also  the  thickness  of 
the  filler  plate,  is  reduced  from  0.0381  m  to  0.0241  m  with  the  introduction  of  bulges  on  both 
sides  to  keep  mass  within  limits.  The  result  shown  in  the  figures  is  for  the  case  Targe  envelope’ 
where  geometric  envelope  limits  are  fairly  large  (relaxed)  on  both  sides.  Lower  envelope  limit 
cannot  be  larger  than  0.3m  due  to  presence  of  the  charge.  With  a  reduced  envelope  limit  of  (-0.1, 
0.14)  m  an  optimum  double  bulge  is  again  seen  but  lies  within  the  envelope.  Table  5  summarizes 
the  results.  The  Targe  envelope’  case  is  slightly  better  compared  to  ‘small  envelope’  case,  as  is  to 
be  expected.  In  Table  5,  both  RMS  displacement,  which  is  the  objective  function  to  be 
minimized,  as  well  as  the  maximum  displacement  are  quoted.  As  was  done  with  the  2-DV  case, 
the  optimized  3-DV  plate  is  also  compared  to  a  flat  plate  of  equal  mass.  Essentially  similar 
reductions  are  again  observed  in  RMS  deflection,  plastic  strain  and  impulse  and  are  omitted  for 
brevity.  Displacement  reaches  its  peak  at  t=  1.1  ms. 
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Plastic  strain  plots  for  baseline  design  and  optimized  design  is  given  in  Figs.  8a,  8b. 
Plastic  strain  is  maximum  at  center  for  the  baseline  design,  as  also  observed  experimentally, 
while  it  is  around  the  borders  of  the  domain  for  the  optimized  panel  (i.e.  it  is  smeared).  This 
demonstrates  better  utilization  of  material.  Maximum  relative  displacement  of  optimized  designs 
versus  baseline  design  is  shown  in  Fig.  9a.  The  total  impulse  responses  of  optimized  designs  are 
compared  to  those  of  baseline  design  in  Fig  9b.  The  saturated  impulse  value  of  Targe  envelope’ 
and  ‘small  envelope’  optimized  panel  designs  are  15.6%  and  13.6%  smaller  compared  to 
baseline  design,  respectively. 


Table  5.  Results  for  baseline  design  and  optimized  designs  for  3-DV  case 


Property 

Baseline 

Design 

Large  Envelope 

(-0.3,1.0)m 

Small  Envelope 

(-0.1,0.14)m 

Optimized 

Design 

Change 

Optimized 

Design 

Change 

Objective  Function  ( m ) 

20.51E-03 

4.34E-03 

78.9% 

4.49E-03 

78.1% 

Max  Relative 

Displacement  ( m ) 

58.43E-03 

7.32E-03 

87.5% 

7.68E-03 

86.9% 

Max  plastic  strain 

0.1277 

0.01932 

84.9% 

0.02028 

84.1% 

Total  Mass  (kg) 

1872.2 

1894.57 

(1.2%) 

1894.54 

(1.2%) 

Saturated  impulse,  (kg- 

m/s) 

6254.3 

5298.62 

15.6% 

5425.82 

13.6% 

15 


Fringe  Levels 
1.277e-01  _ 

1  .ISOe-OI 
1.022e-01 
8.941  e-02 
7.664e-02 
6.387e-02 
5.109e-02 
3.832e-02 
2.555e-02 
1.277e-02 
0.000e+00 

Figure  8a.  Plastic  strain  plot  for  initial  baseline  (flat  plate)  design 


Fringe  Levels 
1.753e-02 
1.578e-02 
1.402e-02 
1.227e-02 
1.052e-02 
8.764e-03 
7.011e-03 
5.258e-03 
3.506e  03 


1.753e-03 

0.000e+00 


Figure  8b.  Plastic  strain  plot  for  3-DV  large  envelope  optimized  design 


16 


Time,  sec 


(a) 


sn 

If 


0> 

v> 

3 

Q. 

E 

rsi 


7000 

ouuu 

5000 

r . 

4000 

snnn 

-♦—Baseline  design 

5UUU 

Large  envelope 

zuuu 

1000 

— a—  small  envelope 

0 


0  2  4  6  8  10  12 

Time,  msec 


(b) 


Figure  9.  (a)  Comparison  of  relative  displacement,  (b)  Comparison  of  impulse  response 


Table  6  compares  the  3-DV  large  envelope  case  to  the  2-DV  large  envelope  case.  While 
the  final  masses  are  a  little  different,  the  3-DV  case  is  better  than  the  2-DV  case  even  with  lesser 
structural  mass. 


Table  6.  Comparison  of  Optimized  Plate  and  Flat  Plate  of  Equal  Mass 


Property 

Flat  plate  of 
equal  mass 

Large  Envelope 
(-0.3,1.0)m 

Small  Envelope 
(-0.1,0.14)m 

Optimized 

Design 

Change 

Optimized 

Design 

Change 

Objective  Function  ( m ) 

1.32E-02 

4.34E-03 

67.12% 

4.49E-03 

65.98% 

Max  Relative  Displacement 

0 m ) 

3.88E-02 

7.32E-03 

81.15% 

7.68E-03 

80.23% 

Max  plastic  strain 

0.08739 

0.01932 

77.89% 

0.02028 

76.79% 

Total  Mass  (kg) 

1894.5 

1894.57 

0.00% 

1894.54 

0.00% 

Saturated  impulse,  ( kg-m/s ) 

6229.3 

5298.62 

14.94% 

5425.82 

12.90% 

Physical  Interpretation  of  Results  of  3-DV  case:  The  bottom  bulge  towards  the  charge, 
observed  in  both  2-DV  and  3-DV  optimized  shapes,  is  likely  playing  an  important  role  in 
deflecting  blast  waves  away  from  the  panel.  The  3-DV  result  achieves  greater  reduction 
compared  to  the  2-DV  result  by  adding  material  to  the  center  while  simultaneously  thinning  the 
plate  to  keep  mass  constant,  what  may  be  termed  a  ‘mass  effect’.  The  amount  of  material  added 
can  be  controlled  by  the  designer  by  suitably  defining  the  envelope  limits. 
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3-DV  Optimization  of  a  plate  supported  on  the  four  edges  on  springs 

Results  for  a  plate  without  the  grip  assembly  and  supported  on  edge  springs  (mimicking 
a  freely  supported  /  lightly  clamped  conditions)  are  presented  in  Table  7.  We  observe,  again,  the 
double  bulge  shape  (Fig.  10)  and  large  improvements  over  the  baseline  design.  However,  as 
noted  in  Section  3,  this  plate  support  system  was  discarded  in  favor  of  the  grip  system  presented 
earlier  owing  to  plastic  strain  concentrations  near  the  corners,  inconsistent  with  experimental 
observations.  The  purpose  of  presenting  this  result  is  to  lend  further  credibility  to  the  optimized 
shaped  obtained  above. 


Table  7.  Optimum  solution  of  plate  with  springs  with  3-DV 


Property 

Baseline 

3DV 

3DV 

(Flat) 

SE 

LE 

Population  Size 

30 

30 

Number  of  Generations 

20 

20 

Max.  relative  displacement  (mm) 

86.2 

31.3 

28.2 

occurs  at  1st  peak 

Objective  =  RMS  displacement, 

46.8 

23.1 

21.8 

mm 

Max.  plastic  strain 

0.055 

0.126 

0.149 

Total  mass  of  structure  (kg) 

152.9 

147.01 

149.87 

Total  Z-momentum  (kN-sec) 

4.93 

3.7 

3.6 

3DV  SE 


3DV  LE 


Figure  10.  Final  optimal  shape  for  plate  with  spring  supports  with  3-DV 
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9-DV  Case:  We  now  present  optimization  results  with  9  basis  shapes.  The  best  shape  obtained 
by  the  optimizer  (repeated  trials  were  needed)  is  shown  in  Fig.  11,  which  is  again  seen  to  be  a 
double  bulge,  albeit  slightly  non-symmetric  in  the  9-DV  case  owing  to  the  optimizer. 
Comparison  with  3-DV  Case  is  given  in  Table  8  for  the  small  envelope  case.  Results  obtained 
show  only  marginal  improvement  over  3-DV  Case.  The  9-DV  Case  is  a  very  challenging  task  for 
any  optimizer.  It  may  still  be  that  a  hitherto  undetermined  unsymmetric  sinusoidal  shape  lies 
within  the  design  space  that  has  a  lower  objective.  However,  we  have  a  fair  amount  of 
confidence  in  the  result  obtained  below  owing  to  the  use  of  various  starting  random  number 
seeds.  Further,  optimal  shapes  that  are  ‘wavy’,  even  if  discovered  in  the  future,  are  hard  to 
manufacture  and  may  be  very  sensitive  to  charge  location. 


Table  8.  Comparison  of  3DV  and  9DV  Cases 


Property 

3  dv 

9  dv 

Objective  function  value  (m) 

4.49E-03 

4.28E-03 

Max  relative  displacement  (m) 

7.68E-03 

7.79E-03 

Max  plastic  strain 

0.020282 

0.02233 

Total  mass  of  the  structure  (kg) 

1894.54 

1895.3 

Impulse  (Ns) 

5401.41 

5313.58 

Figure  11.  Optimum  shape  corresponding  to  9DV  Case 


Sensitivity  to  charge,  standoff  distance:  Limited  sensitivity  studies  were  conducted  to  ensure 
that  the  optimum  shapes  did  not  change  with  ±5  %  and  ±10  %  changes  in  standoff  distance  (z- 
distance  of  charge  from  plate)  and  also  of  charge  amount  (i.e.  kg  of  TNT).  For  brevity,  details 
are  not  given  here,  since  the  optimized  shapes  and  results  are  similar  except  for  small 
perturbations  in  values.  As  the  standoff  distance  decreases,  bulge  towards  the  charge  was  larger 
compared  to  the  amount  of  bulge  on  the  top  face  of  the  panel.  With  the  increases  in  standoff 
distance,  the  optimum  shape  remained  same. 
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8.  Off-Center  Charge  Locations 

The  optimization  is  carried  out  in  four  ‘steps’,  each  successively  offsetting  the  charge  by  a 
distance  of  Lt. 24  (  L  being  the  length  of  the  plate  of  interest)  along  the  x-axis  from  the  initial 
center  position.  The  y,  z  co-ordinates  remain  the  same  at  0,  -0.4064  m  respectively.  In  each 
subsequent  step  the  previous  step’s  optimized  shape  was  used  as  the  initial  shape.  1 1  and  7  basis 
shapes  are  used,  respectively,  for  unsymmetric  and  symmetric  shape  variations.  Shape  morphosis 
as  a  result  of  charge  off-set  is  shown  in  Fig.  12.  Step  3  onwards  a  cavity  appears  on  the  bottom 
surface  of  the  plate  over  the  -x-axis.  The  thickness  (bulge)  of  the  plate  increases  until  the  mass 
limit  is  reached.  Thereafter  a  cavity  is  formed  to  get  more  mass  above  the  charge  location.  The 
C.G.  of  the  plate  shifts  in  the  direction  of  the  off-center  charge.  Figure  13  shows  this  cavity  in 
detail.  Table  9  contains  the  response  values. 


ALUMINUM  PANEL  MINE-BLAST  TEST 


J 


J 


Center  Charge 


Step  1 


Step  2 


J 


J 


Step  3  Step  4 

Figure  12.  Right  Profile  View  of  the  Optimized  Plates  with  increasing  charge  offset 
(towards  the  left,  in  steps  of  L/24  from  center,  L  =  length  of  ‘plate  of  interest’  ) 
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ALUMINUM  PANEL  MINE-BLAST  TEST 

Time  =  0 


Figure  13.  Cavity  appearance  in  Step  4,  at  bottom  of  plate 


Table  9.  Results  with  offset  charge 


Charge  Location 

Objective  =  RMS 
Relative1  Displacement 
(mm) 

Max 

Plastic 

Strain2 

Total  Mass  of 
the  assembly 
(kg) 

Center 

4.48 

0.0231 

1894.54 

Step  1  (L/24) 

4.62 

0.0367 

1895.52 

Step  2  (L/12) 

4.83 

0.0361 

1891.77 

Step  3  (L/8) 

4.91 

0.0293 

1894.25 

Step  4  (L/6) 

5.15 

0.0440 

1893.00 

'Relative  Displacement  =  8  -  8flxtUre  at  1st  peak 
Limit  value  used  is  0.15 

A  verification  study  of  optimality  for  offset  charges:  Owing  to  the  difficulty  in  using  the 
stochastic  optimizer  on  this  problem,  a  coarse  verification  on  optimality  is  carried  out  by  shifting 
the  charge  a  step  backward  (along  x-axis)  and  a  step  forward  for  each  step.  It  is  verified  that  the 
shape  obtained  for  each  step  is  best  for  the  charge  location  for  which  it  was  optimized.  Another 
way  of  interpreting  the  results  is  that  for  a  given  charge  position  the  shape  optimized  for  it  is 
better  than  either  of  the  adjacent  step  shapes.  A  chart  of  the  verification  results  is  provided  below 
in  Table  10. 
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Table  10.  Verification  study  of  optimality  for  offset  charges 
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■  Charge  moved  to 
previous  step 

■  Charge  at  current 
step 

■  Charge  moved  to 
next  step 


9.  Methodology  development  for  optimization  of  honeycomb  sandwich  structure 

Lightweight,  high  strength  and  high  energy  absorbing  materials  are  most  suitable  for 
protective  structures  against  blast  loading.  It  is  well  known  that  metal  sandwich  plates  have 
proved  to  serve  this  purpose  and  honeycomb  sandwich  material  is  a  suitable  option.  For 
maximum  benefit,  different  parameters  of  honeycomb  core  sandwich  metal  plate  need  to  be 
optimized.  However,  modeling  the  each  cell  in  the  honeycomb  core  in  detail  with  iterative 
optimization  techniques  is  computationally  infeasible.  Instead,  the  viable  approach  adopted  here 
is  to  replace  the  honeycomb  core  with  a  solid  plate  with  equivalent  mechanical  properties.  A 
novel  strategy  is  developed  here  to  optimize  the  honeycomb  sandwich  structure  based  on  two 
concepts:  (1)  Parameterization  of  the  unit  cell  in  terms  of  a  finite  number  of  design  parameters,  a 
subset  of  which  are  selected  as  design  variables  for  optimization  after  carrying  out  a  numerical 
study,  and  (2)  Virtual  testing ,  a  novel  technique  developed  here,  to  link  or  characterize  the 
stress-strain  curve  as  a  function  of  the  design  variables,  paving  the  way  for  subsequent 
optimization.  These  two  concepts  are  detailed  below  in  the  context  of  air  blast  mitigation. 

Honeycomb  core  unit  cell  geometry:  The  honeycomb  core  sandwich  structure  considered  here 
consists  of  three  layers:  top  metal  face  plate,  honeycomb  core  (Fig.  14),  and  bottom  metal  face 
plate.  Honeycomb  core  is  characterized  by  its  cell  size  D,  foil  thickness  t,  branch  angel  a,  and 
core  thickness  h.  All  these  parameters  decide  its  crush  strength,  which  is  the  key  property  for 
offering  protection  against  the  blast  load.  In  this  study  we  consider  the  regular  hexagon  cell  core, 
a  =  120°,  as  it  gives  maximum  crush  strength  [15].  The  face  plates  are  characterized  by  their 
thicknesses,  but  these  enter  only  at  the  final  stage  during  optimization  of  entire  structure  and  are 
not  relevant  to  parametrization  of  the  core.  Figure  15  shows  the  triangular  unit  cell  that  has  been 
used  in  this  study.  Figures  15(b)  and  15(c)  show  the  full  model  and  cross  section  of  the  model, 
respectively. 
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Cell  Nfch 


Made 


h  =  Depth 


t  =  Foal  Thickness 


L  =  Longitudinal  Direction 
(Parallel  to  core  ribbons | 

w  =  Transverse  Direction 
( Perpendicular  to  core  ribbons | 


Cell  Sue  |  D) 


Figurel4.  Schematic  view  of  honeycomb  core  geometry 


Figure  15  a  -  c.  Unit  cell  model  of  honeycomb  core 


Virtual  testing  of  Unit  Cell 

A  novel  strategy,  Virtual  Testing,  is  used  to  characterize  the  stress-strain  curve  as  a 
function  of  the  design  variables,  paving  the  way  for  a  subsequent  optimization  study  The  crush 
test  is  carried  out  using  LS-Dyna.  The  foil  is  modeled  by  quadrilateral  Belytschko-Tsay  shell 
elements  and  the  adhesive  layer  of  0.01mm  thick  at  the  double  wall  is  modeled  using  solid  brick 
elements.  A5052  aluminum  alloy  is  used  for  the  foil  in  the  study.  Bilinear  isotropic  material 
model  (Table  11)  is  used  for  the  foil  and  the  adhesive.  Automatic  single  surface  contact  is 
applied  to  the  model  with  sliding  and  sticking  frictional  coefficients  as  0.2  and  0.3  respectively 
[15].  Figure  16  shows  four  different  stages  of  crushing  of  honeycomb  core.  Buckling  of  the  foil 
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(cyclic  folding)  starts  from  near  the  impact  edge  and  propagates  downward  as  the  hammer 
travels.  Figure  17  shows  the  variation  of  nominal  compressive  stress  with  the  volumetric  strain 
(hammer  travel).  Compressive  stress  is  defined  as  the  reaction  force  experienced  by  the  hammer 
divided  by  the  whole  cross  sectional  area  and  the  volumetric  strain  is  calculated  by  the  change  in 
core  depth  divided  by  its  original  value.  The  core  resists  to  buckling  until  the  peak  stress  point 
and  then  undergoes  cyclic  collapse  of  the  foil  as  the  hammer  travels.  The  crush  stress  is  the 
average  of  the  oscillatory  stress  during  the  cyclic  collapse  of  the  foil.  Once  the  entire  core  is 
folded,  then  densification  starts  resulting  in  very  high  compressive  stress.  The  crush  stress  is  a 
vital  property  which  reduces  the  blast  shock  transmission  by  absorbing  the  energy  and  hence 
minimizes  the  damage  to  life.  Although  sufficient  care  is  taken  while  approximating  the  load 
curve,  it  is  not  possible  to  define  the  crush  start  and  end  strain  with  great  accurately. 


Table  11.  Mechanical  properties  of  materials  used  in  the  modeling  of  the  unit  cell 


Material 

Density 

(Kg/m3) 

Young’s 
Modulus  (GPa) 

Yield  Stress 
(MPa) 

Tangent 
Modulus  (GPa) 

Poison’s 

Ratio 

A5052 

2680 

72 

300 

50 

0.34 

Adhesive 

2000 

5 

30 

0 

0.3 

Rigid  metal 

288E5 

200 

- 

- 

0.24 

Figure  16.  Virtual  Testing  simulation  of  gradual  crushing  of  the  honeycomb 
core 
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Volumetric  Strain 


Figure  17.  Load  curve  and  its  different  parameters  of  honeycomb  core 


Effect  of  different  test  parameters  on  load  curve:  Three  different  element  sizes  (Fig.  18)  are 
used  to  study  their  effect  on  the  load  curve.  No  noticeable  change  is  observed  when  element  size 
is  increased  from  0.25mm  to  0.5mm.  Load  curve  doesn’t  shows  any  visible  change  (Fig.  19)  with 
the  change  of  core  depth  from  20mm  to  50mm.  However,  at  a  lower  core  depth  such  as  5mm, 
densification  starts  early.  But  the  crush  stress  is  not  affected  by  the  variation  of  core  depth.  For 
better  design,  the  core  depth  should  be  sufficient  enough  to  allow  few  cyclic  folding  of  the  foil 
during  crushing.  Figures  20  and  21  shows  the  effect  of  the  cell  size  D  and  foil  thickness  t  on  the 
load  curve.  Crush  stress  decreases  with  cell  size  and  increases  with  foil  thickness.  Parametric 
study  is  carried  out  by  varying  the  cell  size  and  the  foil  thickness,  the  element  size  and  the  core 
depth  are  kept  constant  as  0.4mm  and  20mm  respectively. 
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Figure  18:  Effect  of  mesh  size  on  load  curve  Figure  19.  Effect  of  core  depth  on 
load  curve 
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Figure  20.  Effect  of  cell  size  on  load  curve  Figure  21.  Effect  of  foil  thickness  on  load 
curve 


Parameterization  the  load  curve  in  terms  of  design  variables:  Crush  start  and  end  strain 
doesn’t  show  much  change  with  D  and  t,  hence  they  can  be  assumed  to  be  independent.  The  area 
under  the  load  curve  mostly  depends  upon  the  value  of  crush  stress  and  doesn’t  change  much 
with  a  small  error  in  determining  peak  stress,  peak  strain,  crush  start  and  end  strain.  Peak  stress, 
crush  stress  (Fig.  22)  and  final  stress  (Fig.  23)  are  expressed  in  terms  of  dimensionless 
parameter;  t/D  and  trend  lines  are  added.  All  the  curves  fit  well  and  the  relations  match  well  with 
Ref.  [16].  The  tensile  stress  cut  off  is  taken  as  the  average  of  the  final  stress  and  the  crush  stress. 
With  these  relations,  a  load  curve  can  be  expressed  for  a  particular  value  of  D,  t.  Above  analysis 
leads  to  a  load  curve  shown  in  Fig.  24  expressed  in  terms  of  D,  t. 


t/D 


Figure  22.  Parametric  form  of  crush  stress  and  peak  stress 
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t/D 


Figure  23.  Parametric  form  of  final  stress 


Fig-24:  Load  curve  used  for  the  honeycomb  core 


Analysis  of  a  base  line  model  (without  optimization):  Figures  25  and  26  shows  two  types  of 
test  models  of  comparable  weight.  Figure  25  is  the  model  used  earlier  in  the  report,  for  an  all¬ 
aluminum  structure.  No  optimization  is  used  here.  However,  this  analysis  demonstrates  that 
optimization  can  have  significant  benefits.  Material  properties  used  in  the  model  are  shown  in 
Table  12.  Blast  load  of  8  kg  mass  of  TNT  is  applied  from  a  0.4064m  distance  from  the  centre  of 
the  bottom  face  of  the  plate  and  the  results  are  analyzed. 
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Fig-25:  Test  model  for  solid  plate 
sandwich  core 


Fig-26:  Test  model  for  honeycomb 


Table  12.  Material  properties  for  the  model 


Material 

Density 

(Kg/m3) 

Young’s 

Modulus 

(GPa) 

Yield 

Stress 

(MPa) 

Tangent 

Modulus 

(GPa) 

Poison’s 

Ratio 

Tensile 
Stress  Cut 
Off 
(MPa) 

A5083 

2700 

68.9 

225 

0.633 

0.33 

— 

Honeycomb 

Core 

(A5052) 

54.05 

1.28 

0 

6.05 

Steel  4340 

8257.85 

2080 

1010 

0.25 

0.29 

— 

Results  and  discussion  of  analysis:  In  the  sandwich  model,  maximum  z-displacement  and 
effective  plastic  strain  occurs  at  the  bottom  face  plate  (Figs.  27,  28  and  Table  13).  The  blast 
shock  transmits  to  the  honeycomb  core  through  the  bottom  face  plate  and  hence  bottom  face 
plate  undergoes  maximum  displacement  and  plastic  strain.  Most  part  of  the  blast  shock  is 
absorbed  by  the  honeycomb  core,  as  a  result  of  that  top  face  plate  is  protected  (Fig.  29).  The 
basic  idea  of  honeycomb  core  plate  design  is  to  protect  the  structure  facing  towards  the  human 
life.  The  maximum  displacement  of  the  top  face  plate  is  less  than  that  of  plate  model  by  26%. 
The  displacement  of  the  plate  in  plate  model  is  more  oscillating  in  nature  (Fig.  9)  where  as  it  is 
very  smooth  (Fig.  27)  in  the  honeycomb  core  model,  which  indicates  that  major  part  of  the  shock 
is  absorbed  by  the  honeycomb  core.  The  effective  plastic  strain  of  the  top  and  bottom  face  plate 
(Table  13)  is  one  order  less  than  that  of  plate  model.  This  is  because  of  the  fact  that  top  and 
bottom  face  plates  are  loosely  constraint  due  to  easily  deformable  nature  of  the  honeycomb  core. 
This  low  effective  plastic  strain  would  give  us  a  larger  limit  for  shape  optimization  without 
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violating  the  plastic  strain  limit.  Scope  for  size  and  shape  optimization  of  the  honeycomb 
sandwich  structure  is  given  in  Section  10  below. 
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Figure  27.  Max  displacement  of  the  face  plate  Figure  28.  Max  effective  plastic  strain  of 

face  plate 


Table  13.  Result  summary  of  both  the  models 


Model  type 

1st  Peak  relative  Z-displacement 
(mm) 

Maximum  effective  plastic 
strain 

Honeycomb  model 

Bottom  face 

Top  face 

Bottom  face 

Top  face 

114 

43 

0.016 

0.0057 

Plate  model 

58 

0.124 

Figure  29.  Deformation  of  the  honeycomb  sandwich  at  the  1st  peak 
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10.  Conclusions  and  Future  Research 

Shape  optimization  of  a  freely  supported  solid  aluminum  panel  for  air  blast  load 
mitigation  is  carried  out.  Center  and  off-center  charge  locations  are  considered.  This  research  is 
timely  as  very  little  work  has  been  done  in  this  area.  Ls-Dyna  is  coupled  to  a  stochastic  DE 
optimizer  using  a  modular  Fortran  code  that  has  been  developed.  Thus,  accurate  responses  are 
used  at  all  times  during  optimization.  Parallelization  of  the  DE  optimizer  makes  the  optimization 
approach  viable.  The  finite  element  model  has  been  developed  to  reflect  experimental  test 
conditions  and  observed  structural  response.  Two,  three,  and  nine  sinusoidal  basis  shapes  are 
chosen,  respectively.  The  optimum  shape,  a  combination  of  these  basis  shapes,  turns  out  to  be  a 
single-bulge  toward  charge  in  the  2-DV  case  likely  attributable  to  a  deflection  of  blast  waves.  In 
the  3-DV  and  9-DV  cases,  there  is  a  double  bulge,  where  material  is  added  to  the  center  of  the 
plate  with  simultaneous  thinning  of  the  plate.  This  shape  may  be  explained  as  a  combination  of 
deflecting  the  wave  and  a  mass  effect.  The  shape  and  results  are  robust  with  respect  to  small 
changes  in  charge  density  or  standoff  distance.  A  plate  on  springs,  with  no  grip  system  also 
results  in  a  double-bulge  shape. 

Importantly,  at  optimum,  the  plastic  strain  is  smeared,  indicating  better  utilization  of  the 
material.  The  panel’s  RMS  displacement  which  is  the  objective  function,  relative  to  the  fixture, 
is  significantly  decreased  (80%  compared  to  the  baseline  design,  and  67%  compared  to  a  flat 
plate  of  equal  mass).  Saturated  z-impulse  decreased  by  14%.  Maximum  plastic  strain  decreased 
significantly  as  well  (and  smeared  around  the  edges)  and  was  well  within  the  limit.  Results  have 
also  been  verified  using  a  finer  mesh  for  optimized  shapes. 

A  morphosis  of  optimized  shapes  for  increasing  off-set  of  charge  location  show  that  the 
characteristic  double  bulge  shape  remains  with  the  bottom  bulge  shifting  with  the  charge. 
Beyond  a  certain  off-set  value,  the  optimized  shape  changes  in  character,  where  a  bottom  cavity 
is  formed  away  from  the  charge  in  order  to  get  more  mass  above  the  charge  location.  The  C.G.  of 
the  plate  shifts  in  the  direction  of  the  off-center  charge. 

Methodology  for  optimization  of  honeycomb  sandwich  structures  is  developed. 
Geometry  parameters  describing  honeycomb  core  geometry  are  linked  with  the  stress-strain 
curve  and  thereby  to  Ls-Dyna  input  parameters,  using  a  novel  idea  based  on  virtual  testing.  Here, 
virtual  testing  involves  crushing  a  unit  cell  using  Ls-Dyna  again.  Studies  are  conducted  for 
varying  values  of  core  geometry  parameters  leading  to  regression  formulas  that  provide  the 
needed  relationships  between  geometry  and  analysis  inputs. 

The  task  of  finding  a  global  optimum  in  such  highly  nonlinear,  nonconvex  and 
computationally  expensive  functions  is  challenging  and  improved  optimizers  are  needed  in 
future.  The  honeycomb  optimization  methodology  developed  here  needs  to  be  coupled  to 
simultaneous  size  and  shape  optimization  of  the  composite  structure,  involving  thicknesses  of 
face  plates  and  core  depth,  core  geometry,  and  shape  of  outer  boundary  of  face  plates.  The  ability 
of  the  honeycomb  core  to  absorb  energy  reduces  plastic  strain  in  the  top  face  plate  away  from  the 
charge  thus  leaving  significant  room  for  optimization  of  the  sandwich  structure. 

This  work  lays  down  a  methodology  of  structural  shape  optimization  against  blast 

loading. 
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